Riferimenti utili

CRS - coordinate reference system
https://r.geocompx.org/spatial-class#crs-intro
https://r-spatial.org/book/02-Spaces.html#sec-crs

meuse dataset
https://cloud.r-project.org/web/packages/gstat/vignettes/gstat.pdf

metodi di interpolazione geografica
https://r-spatial.org/book/12-Interpolation.html

Obiettivo

Stimare i valori di una variabile nel dominio dello spazio (eventualmente anche del tempo)

Cosa si intende per variabilità spaziale?

• Lo stesso fenomeno ha valori diversi se misurato in punti diversi
• Le differenze quasi mai sono puramente casuali: esiste una dipendenza spaziale (o autocorrelazione)
• I valori dei punti vicini tendono a somigliarsi di più rispetto a punti lontani → questo è ciò che la geostatistica quantifica.

Dati di esempio

Useremo il dataset “meuse” e le librerie che vedete sotto

library(tidyverse) 
library(stats)
library(ggplot2)
library(terra)
library(gstat) 
library(kableExtra)
library(sp)
data(meuse)
class(meuse)

## [1] "data.frame"

kable(head(meuse, 5), caption = "Dataset *meuse*") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Cosa facciamo

Vogliamo un modello per “spiegare” la variabilità e dunque poter stimare (interpolare) al meglio dove non abbiamo informazioni.

Stimare valori di una variabile applicando un modello creato conoscendo valori misurati (e.g. campionamenti). E’ ragionevole adottare un criterio spaziale, ovvero usare un criterio di spazio (il valore del punto sarà probabilmente più simile ad una misura vicina che una misura più lontana, a parità di altre condizioni). Da qui la geostatistica.

fonte: Appunti GIS Prof. Alberto Guarnieri - CIRGEO

Giorno 1

• Richiamo concetti statistici fondamentali: varianza, covarianza, correlazione, autocorrelazione

• Richiamo concetti base su sistemi di riferimento coordinate (CRS), geografiche e proiettate

• Fondamenti statistici dell’analisi spaziale dati e variabili regionalizzate

• Metodi deterministici di stima per interpolazione spaziale (IWD)

Varianza

Vediamo come valutare alcuni momenti di una distribuzione - media e varianza ci danno il primo e secondo momento. Z è la variabile aleatoria, nel caso del dataset MEUSE possiamo scegliere la concentrazione di Zinco.

\[M(Z)=\mu_Z=\mathbb E[Z]\]

\[M(Z)=\mu_Z=\frac{1}{n}\sum_{j=1}^n{Z_j}\]

\[Var(Z)=\sigma^2_Z=\mathbb{E}\Big[\big(Z-\mathbb{E}[Z]\big)^2\Big].\]

\[Var(Z)=\sigma^2_Z=\frac{\sum_{j=1}^n(Z_j-M_Z)^2}{n}\]

Formula più pratica

\[Var(Z)=\sigma^2_Z=\mathbb{E}[Z^2]-\mathbb{E}[Z]^2\]

\[Var(Z)=\sigma^2_Z=M(Z^2)-M(Z)^2 \]

Varianza esempio

Calcoliamo la varianza con tre formule diverse:

MZ<-mean(meuse$zinc) 
VarZ <- var(meuse$zinc)
## media quadratica dei quadrati
MZ2<-mean(meuse$zinc^2) 
VarZ <- MZ2 - MZ^2
VarZ <- mean((meuse$zinc-MZ)^2)
  
 

plot( meuse$y, meuse$zinc, xlab="Coord Y", ylab="Zinco", pch="+")
abline(a = MZ,0 )
rect(min(meuse$y), MZ-VarZ^0.5, max(meuse$y), 
     MZ+VarZ^0.5, col = "#ff000045", border = "#ffffff00")

# varianza  = (dev std)^2

Varianza esempio (2)

hist(meuse$zinc)
rug(meuse$zinc)
rect(MZ-VarZ^0.5, 10000, 
     MZ+VarZ^0.5, 0, col = "#ff000045", border = "#ffffff00")
abline(v=MZ, col="red", lwd=2)

Covarianza

La covarianza tra due variabili casuali \(Z_1\) e \(Z_2\) misura quanto esse tendono a variare insieme. Se positiva variano entrambe nella stessa direzione, se negativa in direzioni opposte.

\[\mathrm{Cov}(Z1,Z2)=\mathbb{E}\Big[\big(Z1-\mathbb{E}[Z1]\big)(Z2-\mathbb{E}[Z2]\big)\Big]\]

Esempio Covarianza

La covarianza dipende dalle unità di misura delle variabili → non è standardizzata. Per questo si usa spesso la correlazione, che normalizza la covarianza tra -1 e 1.

La correlazione = 1 equivale ad una covarianza che ha valore uguale alla varianza.

CovZ<-cov(meuse$zinc, meuse$zinc)

var(meuse$zinc) == CovZ

## [1] TRUE

Covarianza

Analizziamo un gradiente di temperatura simulato

load("oggetti.rda")
plot(lat_seq, temps, type = "l", col = "blue", lwd = 2,
     xlab = "Latitude (°N)", ylab = "Temperature (°C)",
     main = "Gradiente Temperatura: Italia → Danimarca")
points(lat_seq, temps, pch = 16, col = rgb(0,0,1,0.3))

Covarianza

La covarianza equivale alla varianza se la calcoliamo confrontandola a se stessa, ma se trasliamo di un dato intervallo (“lag”) e la ricalcoliamo, la covarianza non diventa zero, diminuisce gradualmente….

cov(temps,temps)

## [1] 6.150969

cov(temps[2:200],temps[1:199])

## [1] 5.122209

cov(temps[3:200],temps[1:198])

## [1] 4.972588

## Provate con una variabile senza correlazione 
# rand <- runif(200) 
# cov(rand,rand)
# cov(rand[2:200],rand[1:199])
# cov(rand[3:200],rand[1:198])

Questo è un concetto molto importante da comprendere.

Covarianza

Più aumentiamo la distanza del lag, come è ragionevole pensare per dati spazialmente correlati, più diminuisce il valore di covarianza (si allontana dal valore della varianza == covarianza ).

Se si abbassa per poi rialzarsi, può indicare un pattern ricorrente (provate sostituendo i valori di temps con un valore ciclico come temps <- sin((1:200)/10) .

nm <- length(temps)
plot(temps[1:200], temps[1:200], pch="x", cex=0.5, xlab="Temperatura", ylab="Temperatura" )
points(temps[1:199], temps[2:200], pch="x", cex=0.5, col="red" )
points(temps[1:198], temps[3:200], pch="x", cex=0.5, col="green" )

Covarianza

Non ci crediamo? Rifacciamolo per un dato casuale con media e varianza uguale a quella del nostro dato:

t.rand <- rnorm(nm, mean = MZ, sd = VarZ^0.5)

plot(t.rand[1:(nm-1)], t.rand[2:nm], pch="x", cex=0.5, xlab="km", ylab="Temperatura" )
cc=1
covs.rand<-list("0"=cov(t.rand, t.rand), "1"=cov(t.rand[1:(nm-1)], t.rand[2:nm] ) )

for(lag in seq(6, 36, 10)){ 
  covs.rand[[as.character(lag)]] <- cov(t.rand[1:(nm-lag)], t.rand[(lag+1):nm])
  points(t.rand[1:(nm-lag)], t.rand[(lag+1):nm], cex=0.5, pch="x", col=cc )
  cc=cc+1
}

plot(names(covs.rand),  unlist(covs.rand), type="b", 
     ylab="Covarianza", 
     xlab="lag" )

Correlazione

Relazione lineare tra due variabili, cosiddetta di Pearson. Semplicemente normalizza la covarianza per le deviazioni standard delle due variabili moltiplicate tra loro:

\[-1\le\rho_{Z1Z2}=\frac{Cov(Z1Z2)}{\sqrt{Var(Z1)}\cdot\sqrt{Var(Z2)}}=\frac{\sum_{i=1}^n(Z1-{\mu_{Z1}})(Z2-{\mu_{Z2}})}{\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z1}})^2}\sqrt{\sum_{i=1}^n(y_i-{\mu_{Z2}})^2}}\le+1\]

Autocorrelazione

L’autocorrelazione misura quanto una variabile è correlata con sé stessa quando viene considerato uno spostamento nel dominio dei valori.

In una sola dimensione — sia essa tempo o spazio (ad esempio nel nostro transetto di temperatura SUD → NORD) — questo spostamento viene detto lag, traducibile come sfasamento.

L’analisi dell’autocorrelazione si effettua tramite la funzione CCF (Cross-Correlation Function), che nel caso specifico valuta la correlazione della serie con sé stessa a diversi lag.

acf() in R ti mostra l’autocorrelazione ai vari lag (sfasamenti).

• Lag = 1 → confronta ogni punto con il successivo.

• Lag più grandi → confrontano punti sempre più distanti lungo il transetto.

#ccf(temps, temps, plot = F)
#ccf(temps, temps)

# Autocorrelazione
acf(temps, main = "Autocorrelazione della temperatura")

Cosa sono i dati spaziali?

Rappresentare informazioni associate alla loro posizione nello spazio

Due modelli: il modello vettoriale ed il modello raster.

Ogni valore è georiferito, ovvero associato ad un punto in un sistema di riferimento, quindi delle coordinate.

L’accuratezza può variare a seconda dello strumento usato per la misura.

L’accuratezza dipende dalla scala, che determina la tolleranza agli errori.

Modello vettoriale

I dati vettoriali hanno tre primitive geometriche di base: punti linee e poligoni. La base è sempre il punto con una coordinate, in 2d o 3d.

By M. W. Toews - Own work, CC BY-SA 4.0

Ci sono decine di formati digitali per i dati vettoriali.

Per un riferimento vedete la pagina QUI.

Ad oggi gli strumenti in ambiente R di riferimento per leggere dati vettoriali è la libreria sf.

Modello raster

Il formato raster è una semplificazione/discretizzazione dello spazio in celle/pixel/ nodi a passo costante, quindi una “griglia” georiferita mediante uno o più punti. Il passo della griglia è spesso identico nelle due direzioni (ma non necessariamente), e ne determina la risoluzione spaziale.

fonte: appunti GIS Prof. Alberto Guarnieri

Modello raster

Può essere 3D (voxel)

Ci sono decine di formati di file che rappresentano dati raster. Per un riferimento vedete la pagina QUI.

Ad oggi la libreria di riferimento per leggere questi dati è la libreria terra, e stars

Attenzione alla gestione della memoria in quanto i dati raster possono diventare molto pesanti dal punto di vista di spazio di memoria (1000 pixel x 1000 pixel è un raster “piccolo” ma già fornisce 1e6 elementi.

Sistemi di riferimento

Sono delle convenzioni - si stabilisce un sistema con una origine ed un orientamento degli assi tra loro perpendicolari - Coordinate planimetriche o polari - Sistemi locali o globali

Sistemi di riferimento - geoide

Prima di tutto bisogna definire un riferimento fisico - la forma della superficie definita come “interfaccia solida” è “complicata” (è elastica e subisce l’effetto di svariate forze che la deformano)… e neanche stabile nel tempo!

Si prende come convenzione la forma che avrebbe la superficie equipotenziale (livello del mare) se non ci fossero terre emerse - lo si chiama GEOIDE. E’ la superficie di riferimento per comprendere dove la gravità spinge l’acqua. NB è in corso con il recente piano di telerilevamento nazionale il calcolo di un nuovo nuovo geoide ufficiale nazionale (vedi qui)

Sistemi di riferimento - ellissoide

Il geoide è fondamentale per comprendere fenomeni legati alle quote, ma molto complesso da trattare. Si semplifica con l’ellissoide, una rappresentazione trattabile matematicamente.

Ci sono tanti ellissoidi che sono stati stimati per avvicinarsi il più possibile alla realtà del geoide

Coordinate geografiche

Latitudine e longitudine, coordinate polari.

Coodinate geocentriche

Coordinate cartografiche

Coordinate planimetriche, su un piano. Ma la terra non è piatta!

Soluzione: accettiamo deformazioni ma le teniamo sotto controllo (scala, accuratezza, tolleranza)

Coordinate cartografiche

Si proiettano i punti nella “sfera” su un piano.

Fonte: QGIS documentation

I sistemi di riferimento (CRS)

Vedi documentazione di QGIS QUI

Vedi NOTA PER IL CORRETTO UTILIZZO DEI SISTEMI GEODETICI

Un esercizio con R

Le coordinate di Prato della Valle: Lat: 45.398557, Long. 11.876460 trasformate nei vari sistemi. TEST - prova a consultare il documento IGM e trasformare la coordinata nel sistema “RDN2008 / Zone 12 (N-E)” e vedere come differisce dal sistema “RDN2008 / Italy zone (N-E)” … sono uguali?

library(sf)
lat<-45.398557
long<-11.876460
# punti in classe SF
pt <- st_sfc(st_point(c(long, lat)), crs = 4326)
# UTM fuso 32
pt_utm <- st_transform(pt, 32632)
# 
pt_ecef <- st_transform(pt, 4978)
# le coordinate nei vari sistemi
geographic = sf::st_coordinates(pt)
planimetric = st_coordinates(pt_utm)
geocentric = st_coordinates(pt_ecef)

Variabili regionalizzate

variabili misurate in punti e correlate spazialmente

• Possono essere trattate come campioni di un fenomeno continuo (non discreto)

• Presentano autocorrelazione spaziale

Statistica descrittiva di base

Vedi slides iniziali

• Media, mediana, varianza, deviazione standard

• Distribuzione dei dati (istogrammi, boxplot)

• Identificazione di outlier - per il trattamento rigoroso delle osservazioni (non lo trattiamo)

• Differenze tra dati indipendenti e dati spazialmente correlati

Metodi deterministici vs stocastici

I metodi deterministici non considerano la componente aleatoria (incerta, casuale) dei dati, ma si basano su assunzioni deterministiche per stimare i valori in punti non campionati - ovvero non considerano la struttura spaziale, nessuna autocorrelazione e non forniscono incertezza.

Stimano il valore Z* in un punto dove non abbiamo una misura usando le misure note vicine, ma come?

\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso del punto \(1\) dipende dalla distanza \(h\) di un \(i-esimo\) numero di punti \(n\)

Vicino più prossimo (Nearest Neighbour - NN)

Semplicemente quello più vicino
\[ Z*= Z_i \] Dove \(Z_i\) è il vicino più prossimo e \(Z_*\) il punto da interpolare

Vicino Naturale (Natural Neighbour - NatN)

poligoni di Thiessen danno l’area di influenza.

\[ Z*= Z_i \] dove il peso è indicato da Sibson

\[ Z_i(\mathbf{x})=\frac{A(\mathbf{x}_i)}{A(\mathbf{x})} \]

dove \(A(x)\) è il volume della nuova cella centrata in x e A(xi) è il volume dell’intersezione tra la nuova cella centrata in x e la vecchia cella centrata in xi.

By Markluffel - Own work, CC BY-SA 3.0

• interpolatore esatto: i valori dei dati originali vengono mantenuti nei punti di riferimento.

• crea una superficie liscia priva di discontinuità.

• è interamente locale, in quanto si basa su un sottoinsieme minimo di posizioni dei dati che esclude le posizioni che, pur essendo vicine, sono più distanti di un’altra posizione in una direzione simile.

• è adattivo dal punto di vista spaziale e si adatta automaticamente alle variazioni locali della densità dei dati o della disposizione spaziale.

• privo di ipotesi statistiche (puramente deterministico).

• può essere applicato a insiemi di dati molto piccoli, poiché non è basato su dati statistici.

• privo di parametri, quindi semplice

Inverso pesato sulla distanza (IWD)

\[Z* = F(x,y) = \sum_{i=2}^n w_i Z_i\] dove il peso dipende dalla distanza \(h\)

\[w_i = \frac {h_i^{-p}}{\displaystyle \sum_{j=i}^n h_j^{-p}}\]

e da una potenza \(p\) che va inserita dall’utente e può essere oggetto di ottimizzazione.

By Hassan Hadji CC-BY-SA Wiki

Esempio IDW

library(gstat)
data(meuse.grid)
data(meuse)
## molto importanti i passaggi sotto
coordinates(meuse) = ~x+y
coordinates(meuse.grid) = ~x+y
gridded(meuse.grid) = TRUE
### 
zinc.idw = idw(zinc~1, meuse, meuse.grid)

## [inverse distance weighted interpolation]

class(zinc.idw)

## [1] "SpatialPixelsDataFrame"
## attr(,"package")
## [1] "sp"

spplot(zinc.idw["var1.pred"], main = "Interpolazione IDW concentrazione Zinco")

Esportiamo su raster GeoTIFF per GIS

La libreria terra riconosce l’oggetto “sp” - posso aprirlo con qualsiasi sofware GIS! Nota: non vengono specificati sistemi di riferimento.

zinc.idw.terra <- terra::rast(zinc.idw)
plot(zinc.idw.terra, main="Mappa Concentrazione Zinco")

# terra::writeRaster(zinc.idw.terra, "zincIDW.tif")

Esportiamo su raster GeoTIFF per GIS

Conoscendo l’area, vicino a Maastricht - Paesi Bassi - applichiamo un SR probabile e verifichiamo aprendo in GIS

terra::crs(zinc.idw.terra) <- "EPSG:28992"
# terra::writeRaster(zinc.idw.terra, "zincIDW.tif", overwrite=T)

Esercizio finale giorno 1

Eseguite la procedura sopra per interpolare la variabile “lead” (piombo). Salvate il GeoTIFF con il sistema di riferimento corretto nella cartella condivisa QUI dentro la cartella “giorno1” nominando il file cognome_nome.tif